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Abstract The spectra of low-lying pair excitations for an imbalanced two-component 
superfluid Fermi gas are analytically derived within the path-integral formalism 
taking into account Gaussian fluctuations about the saddle point. The spectra are 
obtained for nonzero temperatures, both with and without imbalance, and for ar- 
bitrary interaction strength. On the basis of the pair excitation spectrum, we have 
calculated the thermodynamic parameters of state of cold fermions and the first 
and second sound velocities. The parameters of pair excitations show a remark- 
able agreement with the Monte Carlo data and with experiment. 
PACS numbers: 03.75.Ss, 05.30.Fk, 03.75.Lm 



1 Introduction 

Recent experimental breakthroughs in the manipulation of ultracold Bose and 
Fermi gases have opened new prospects for advancing many-body physics 1 . The 
dimensionality of these gases can be controlled with optical lattices, and the inter- 
action strength can be tuned using Feshbach resonances. The experimental control 
over geometry and interactions in ultracold atomic gases has turned these systems 
into powerful quantum simulators that can test and generalize many-body theories 
originally developed for solid state systems. Recently, much attention has been 
paid to ultracold atomic gases with strong interactions because of their possible 
relation to some striking natural phenomena including high-temperature super- 
conductors and neutron stars^. 

In particular, great experimental and theoretical effort has been devoted to the 
study of superfluidity arising from pairing in ultracold Fermi gases 3 . Specifically, 
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the effect of population imbalance (between the pairing partners) on the super- 
fluid pairing mechanism is a topic of current investigation. Of no less interest in 
the context of high-temperature superconductivity is the study of the crossover 
between the Bardeen-Cooper-Schrieffer (BCS) and the molecular Bose-Einstein 
condensation (BEC) regimes. 

In order to describe the superfluid phase transition as well as the broken- 
symmetry phase below a critical temperature T c , several methods based on the 
T-matrix approach have been develope d 4 ' 5 ! 6 ' 7 ' 8 ' 9 ! 10 . Amongst those, the Nozieres- 
Schmitt-Rink (NSR) theory 4 and its path-integral reformulation 5 have been very 
successful and remain widely used. Here, we use an improved version of the NSR 
scheme, that we will denote as the Gaussian pair fluctuation theory (GPF) U ' 12 , 
and that effectively works both at low temperatures and above T c . 

As shown in Ref. 13 , the thermodynamic properties of the superfluid Fermi gas 
at sufficiently low temperatures can be derived within a simple model using the 
spectrum of low-lying elementary excitations. In that model, the thermodynamics 
of the superfluid Fermi gas 13 is based on the fermion-boson model with phe- 
nomenological parameters whose values are determined from the Monte Carlo 
calculation s 141 °i 18 i 19 (i n the limit of zero temperature). The model exploited in 
Ref. 13 describes well the thermodynamic properties of a balanced unitary Fermi 
gas at low temperatures. Moreover it was demonstrated 20 that in the BEC limit, 
the imbalanced Fermi superfluid indeed reduces to a simple Bose-Fermi mixture 
of Bose-condensed molecules and unpaired fermions. The goal of the present 
work is to extend these results to non-zero temperatures, and to arbitrary scat- 
tering lengths. 

For this purpose, we calculate the parameters needed for the fermion-boson 
model from the NSR and/or GPF theories. Our results are compared with the 
Monte Carlo data at unitarity. We analytically derive the spectra of low-lying el- 
ementary excitations for an imbalanced Fermi gas in 3D at finite temperatures in 
the whole range of the BCS-BEC crossover. These spectra are obtained using the 
path-integral representation— of the NSR theory extended to imbalanced Fermi 
gase s 21 ' 22 as well as the GPF approac h n i 12 . Using the obtained spectra of the el- 
ementary excitations, thermodynamic parameters such as the internal energy, the 
chemical potential, the first and second sound velocities, are calculated. 

2 Formalism 

2.1 Path-integral GPF approach for imbalanced Fermi gases 

We consider a two-component Fermi gas within the path-integral approach. The 
path-integral formulation 5 of the NSR scheme has been extended in Refs. 21 ' 22 
to the case of unequal 'spin up' and 'spin down' populations of fermions. In the 
present work, the treatment of the imbalanced Fermi gas is performed using the 
NSR scheme 21 ' 22 and its improved version, the GPF theory 11 extended to the 
imbalanced case. 

The thermodynamic parameters of the imbalanced Fermi gas are completely 
determined by the thermodynamic potential Q of the grand-canonical ensemble. 
The thermodynamic potential Q, the same as in Refs. 21 ' 22 , is the sum of the 
saddle-point thermodynamic potential Q sp and the fluctuation contribution Qp. 
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These thermodynamic potentials are provided, respectively, by the zeroth-order 
and quadratic terms of the expansion of the Hubbard-Stratonovich pair-field ac- 
tion around the saddle point. 

The saddle-point thermodynamic potential for the imbalanced Fermi gas with 
i-wave pairing is- 2 -^ 



sp 
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where V is the system volume, /3 is the inverse to the thermal energy k^T, A is 
the amplitude of the gap parameter, a s is the scattering length, ^ = k 2 — /J, is the 

fermion energy, and = y^+A 2 is the Bogoliubov excitation energy. The 

chemical potentials of imbalanced fermions are expressed through the averaged 
chemical potential fi = (jj,^ + jJL^j /2 and the chemical potential imbalance £ = 
(jUf — jU|) /2. We choose the units with h = 1, the fermion mass m = 1/2, and 

the Fermi energy Ef = h 2 (3n 2 n) 2 ^ / (2m) = 1 (n is the total fermion density). 
The fluctuation contribution to the thermodynamic potential is the same as in 
Refs.S 

The gap parameter is found from the gap equation minimizing the saddle -point 
thermodynamic potential, 
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from which we can extract A(T,/J,,Q. For an imbalanced gas, the saddle-point 
thermodynamic potential can have two minima: one at A = and one at A ^ 0. 
This can result in a first-order superfluid phase transition 23 . With our notation, we 
emphasize that the thermodynamic potential is a function of T, jJ. , £ (and actually 
V, but this dependency drops out). However, A is treated as an additional parame- 
ter on which the thermodynamic potential depends. It is this treatment of A as an 
additional parameter (in the broken-symmetry phase with 4^0) that leads to a 
distinction between the NSR approach and the GPF approach. When calculating 
the gap equation one should use 
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In the standard NSR approach the last terms (involving the derivatives of A) are 
omitted. The GPF approach suggested in Refs. 1112 takes into account the addi- 
tional derivatives for the balanced case, and corrects the NSR densities for changes 
in A as fi and £ are varied. The GPF method presented here is the path-integral 
formulation of the GPF theory of Refs. n i 12 extended to the imbalanced case. As 
shown in Ref. 12 , the GPF theory provides the best overall agreement of its ana- 
lytic results with experiment and with Monte Carlo data, except in close vicinity 
to T r . 
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The GPF corrections are not present when one uses the saddle-point approxi- 
mation and calculates 



l sp 



8n sp = - 



dQ, p (T,fi,Z;A) 



(5) 
(6) 



However, the correction terms will be relevant for the calculation of the fluctuation 
contributions, nf\=n — n sp and 8iifi = 8n — 8n sp . These fluctuation contributions 
to the density and Srifi are given by the same expressions as in Ref. 21 
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Here, no is an arbitrary positive integer, and the parameter y lies between two 
bosonic Matsubara frequencies V„ < y < V„ 0+ i, V„ = 2nn/j3. The spectral func- 
tions J (q,z) and (q,z) of complex frequency z are 



J(q,z) = 
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where M;* (q,z) are the matrix elements of the pair field propagator. The matrix 
elements M;* (q,z) are given by the expressions 21 
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Here, the distribution function X (E k ) is 
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The derivatives in (0 and d 1 Ob within the GPF scheme are determined as men- 
tioned above - taking into account a variation of the gap parameter. Within the 
NSR scheme, these derivatives are calculated assuming A to be an independent 
variational parameter. 

The equation of state of the imbalanced Fermi gas is thus determined as a joint 
solution of the saddle-point gap equation and the number equations accounting 
for Gaussian fluctuations. Within both the NSR and GPF schemes, the Gaussian 
fluctuations do not feed back into the saddle -point gap equation. 



2.2 Low-lying pair excitations 

In order to obtain the spectrum of low-lying pair excitations for the imbalanced 
Fermi gas, we perform the long-wavelength and low-energy expansion of the ma- 
trix elements M;j (q,z) as proposed in Ref. 10 . We take into account the terms up 
to quadratic order in powers of q and z, and find: 



The coefficients of the expansion (114b are derived straightforwardly. After 
some algebra, we arrive at their expression through the integrals: 
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The low-lying pair excitations correspond to the poles of the spectral func- 
tions (O and d 1 Ob . Therefore the dispersion equation for the energies of the pair 
excitations ft) = (0 q is 

M M (q, ffl)Mi,i (q, -0)) -M 2 2 (q, ft)) = 0. (16) 

This equation is solved expanding ft) 2 up to the terms of the order of q 4 . We then 
obtain the energies C0 q in a form structurally similar to the collective excitations in 
Ref-i 3 -: 

(0 q = y^V + V, (17) 

where the parameters c and A are related to the coefficient of the expansion (1141) 
as follows, 

/ 2A(B-E) \ 1/2 
C -U 2 + 2A(//-G)J ' (18) 
x = C 2 (B- E) (44 (gg - Eg) + C 2 (g + £)) 
(C 2 +2A(//-g)) 3 

For small pair momentum q, the energy of the pair excitation becomes linear 
in the momentum. Thus (O q at small q represents a Bogoliubov- Anderson mode, 
which is gapless in accordance with the Nambu-Goldstone theorem. The param- 
eter c has the dimensionality of velocity. In the zero-temperature limit for a bal- 
anced gas, all fermions are in the superfluid state, and the velocity parameter for 
the pair excitations tends to the first sound velocity for the whole fermion system. 
The so-called gradient parameter A provides a growth of kinetic energy due to a 
spatial variation of the density 24 ' 25 . 

The pair excitation spectra obtained in the present work generalize the long- 
wavelength expansion of Ref. 10 to the case of non-zero temperatures and unequal 
'spin-up' and 'spin-down' fermion populations. 

In Fig. 1, the parameters c and A characterizing pair excitations are plotted 
as a function of temperature for the balanced Fermi gas in the unitarity regime 
(l/a s =0). The parameters of the equation of state for a given temperature are 
determined from the joint solution of the gap and number equations taking into 
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Fig. 1 (color online) (a) Sound velocity parameter of the Bogoliubov-Anderson mode c for cold 
fermions in 3D as a function of temperature in the unitarity regime and for £ = 0. The parameters 
of state are determined taking into account the fluctuations in different ways: within the path- 
integral GPF formalism (solid curve) and within the NSR scheme (dashed curve). The symbols 
show the results of the different Monte Carlo calculations 14 i 15 i 16 i 17 [ 18 . (ft) The parameter X as 
a function of temperature in the same regime. The dot-dashed vertical line indicates the critical 
temperature of the superfluid phase transition. The thin dotted curves show the formal solutions 
for the parameter c and X within the path-integral GPF and NSR approaches above T c . Inset: 
Chemical potential calculated within the path-integral GPF and NSR approaches compared with 
the Monte Carlo data from Ref 



account fluctuations in the number equation. The fluctuation contributions to the 
fermion density are calculated within the path-integral GPF and NSR schemes. 

As found in Ref. 1 1 , the NSR approach becomes inaccurate in the vicinity of 
the critical temperature T c . It was also shown that the NSR scheme reveals a re- 
entrant behavior of the parameters in the state above T c , leading to an artificial 
first-order superfluid phase transition—. The re-entrant behavior of the parame- 
ters c and A obtained in the NSR approach is clearly seen in Fig. 1 . The critical 
temperature T c for the balanced gas, indicated by a dash-dotted line in Fig. 1, is 
the same within the NSR and GPF approaches. However, the GPF method leads 
to better results with respect to NSR for the broken-symmetry phase. This can be 
seen, for example, in the inset of Fig. 1 where we plot the chemical potential as 
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a function of temperature below T c calculated within the NSR and path-integral 
GPF approaches and compared with the Monte Carlo results of Ref.—. 

In the zero-temperature limit, the sound velocity parameter c for the pair exci- 
tations obtained within both the path-integral GPF and NSR approaches exhibits 
an excellent agreement with the numerical results obtained using different Monte 
Carlo algorithm s 14 i 15 i 16 i 17 i 18 , Also the gradient parameter A at zero temperature 
lies within the range of the values of A obtained in Refs. 13 ' 25 as the best fitting 
parameters for the ground state energy of fermions compared with Monte Carlo 
data. This agreement demonstrates the accuracy of the present approach for the 
broken-symmetry phase of cold fermions. 
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Fig. 2 (color online) Sound velocity parameter (a, c) and parameter X (b, d) for the Bogoliubov- 
Anderson mode as a function of the inverse scattering length for kgT = O.OlEp (solid curves), 
for kgT = 0.05Ef (dashed curves) and for kgT = O.lEf (dotted curves) for a balanced Fermi 
gas (a, b) and at the chemical potential imbalance £ = 0.2 (c, d) in the unitarity regime. 



In Fig. 2, we plot the parameters c and k as a function of the inverse scattering 
length 1 Ja s for the balanced gas (the left-hand panels) and at the chemical poten- 
tial imbalance £ = 0.2 (the right-hand panels). At kgT = O.OlEf the sound veloc- 
ity parameter monotonously decreases with increasing \/a s . For non-zero temper- 
atures, c exhibits a maximum, which shifts to higher coupling strengths for higher 
temperatures. The parameter X at finite temperatures has a minimum, which al- 
most vanishes in the zero-temperature limit. In the weak-coupling regime, both 
c and A are sensitive to temperature. When moving towards the strong-coupling 
regime, c and A gradually become almost independent on T. The imbalance leads 
to the appearance of a critical inverse scattering length such that for smaller 1 /a s , 
there is no superfluid state (see, e. g., Ref.—). 
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3 Parameters of state 

3.1 Thermodynamic functions 

Using the spectra of the elementary fermionic and pair excitations derived in Sec. 
[21 we can obtain the thermodynamic functions of the superfiuid Fermi gas at finite 
temperature. In Ref. 13 , a similar description of the thermodynamic properties was 
performed for a unitary balanced Fermi gas using the zero-temperature spectra of 
elementary excitations. 

In Ref. 13 , a pair excitation spectrum of the form of expression ( 117b is used, 
where the zero-temperature sound velocity c is taken from the Monte Carlo dat a 14 i 15 ' 16 i 17 > 18 
and the gradient parameter X is determined from a fit of the thermodynamic prop- 
erties to the Monte Carlo results. In the present calculation, the pair excitation 
spectra are obtained using the analytic path-integral GPF approach without any 
fit. 

In this section we consider the thermodynamic functions of the cold Fermi gas 
within the model of fermionic and pair excitations. The grand-canonical thermo- 
dynamic potential is the sum of the saddle-point and pair excitation contributions 

Q =Q sp + Q p . (20) 

The saddle-point thermodynamic potential Q, sp is given by Eq. (QJ. The contribu- 
tion of pair excitations isi^ 

The entropy S is found using its relation to the grand-canonical thermody- 
namic potential, 

Using the thermodynamic potential Q with (Q]) and (I2TI) we find that the entropy 
is expressed as 

s = Vjf_ r & E k (cosh/3E,cosh/3g + 1) - gsinh^sinhffg ^ 



3?r 2 Jo E k (coshjSC + cosh/^) 2 

V r°° I g~P m q / \\ 

+ 2^L (^izt^-K 1 -«~'*)J ^- (23) 

Finally, the internal energy E of cold fermions is calculated using the relation 
between E and the grand-canonical thermodynamic potential, E = Q + TS + 
jj.N, where is the total number of fermions. 

In Fig. 3, the internal energy per particle for a unitary Fermi gas calculated 
in different approaches is plotted as a function of temperature for the broken- 
symmetry phase at T $J T c . The critical temperature determined within the path- 
integral GPF model is the same as within NSR, T c sa 0.225Ep /kg. The results of 
the present calculation within the model of fermionic and low-lying pair excita- 
tions with parameters determined using the path-integral GPF and NSR methods 
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are shown with short-dashed and dot-dashed curves, respectively. The other re- 
sults represented in Fig. 3 are (after Ref. 13 ): the internal energy calculated within 
the low-temperature fermion-boson (FB) model 13 , and the result of the analytic 
model proposed by Bulgac, Drut, and Magierski (BDM) 19 . The analytic results 
are compared with those of the Monte Carlo calculations from Ref. 19 and with the 
experimental data of Ref. 27 for a gas of 6 Li atoms at unitarity. 
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Fig. 3 (color online) Internal energy per particle for a Fermi gas at unitarity calculated within the 
model of fermionic and low-lying pair excitations with parameters determined using the path- 
integral GPF method (short-dashed curve) and the NSR theory (dot-dashed curve). The solid 
curve shows the internal energy obtained within the low-temperature fermion-boson model in 
Ref—. The dashed curve is the result of the BDM model 19 . The full dots represent the Monte 
Carlo simulations 19 . The squares are the experimental data—. 



As reported in Ref. , the low-temperature fermion-boson model works well 
in the broken-symmetry phase where the internal energy resulting from this model 
is close to the Monte Carlo data of Ref. 19 . The present study is performed with the 
parameters of the elementary excitations obtained using the analytic approaches 
rather than a fit to the Monte Carlo simulations. The internal energy calculated 
in the present work with the parameters determined using the path-integral GPF 
method is very close to the Monte Carlo results at low temperatures T J$ 0.55T C . 
Furthermore, our result is in good agreement with the experiment 27 in the whole 
temperature range below T c . 

Neither the model of fermionic and pair excitations used in the present work, 
nor the low-temperature fermion-boson model of Ref. 13 predicts the superfiuid 
phase transition: the critical temperature T c is determined within the path-integral 
GPF method before the long-wavelength expansion is performed in Sec.|2j How- 
ever, the latter model can describe well the broken-symmetry phase of cold fermionic 
atoms. 
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3.2 Superfluid density 

The total density of cold fermions within the model of fermionic and pair excita- 
tions is given by the sum of fermion and boson contributions 

n = i j>* ( l - k Xk ) + (24) 

The total density is a sum of the normal and superfluid densities n = n n + n s . The 
superfluid density n s , as well as the total density, is constituted by the saddle-point 
result for an imbalanced Fermi gas and the contribution due to the pair excitations, 

+ ir(^-f^ 2 p^)^ (25) 

where the function is given by 

_ dXk cosh cosh /3£ + 1 
_ dE k * 

(cosh jS£7yt + cosh j3 C, )" 

(see the analogous expression for the Fermi gas in 2D in Ref. 28 ). In the balanced 
case, the superfluid density becomes equivalent to the corresponding expression of 
Ref. 13 . but with other values of the parameters of the pair excitations, as discussed 
above and shown in Fig. 1 . 

In Fig. 4, the superfluid density divided by the total density calculated within 
the model of fermionic and pair excitations is plotted as a function of temperature. 
As seen from the figure, the superfluid density calculated using the parameters of 
the pair excitations obtained within the NSR model exhibits a re-entrant behavior 
above T c similarly to the sound velocity parameter in Fig. 1 . The analogous bend- 
over of the superfluid density above T c was reported in the full NSR approach in 
RefJS. 



3.3 Sound velocities 



We consider the sound propagation in a superfluid Fermi gas using the approach 
of the two-fluid hydrodynamic s 29 ' 30 in the same way as in Ref. 13 . The first sound 
velocity u\ in the two-fluid hydrodynamics is determined by the formula 
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(26) 
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with the entropy per particle S = S/N. The pressure is proportional to the grand- 
canonical thermodynamic potential: P = Q /V. We adopt the expression— 
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Fig. 4 (color online) Superfluid density for a Fermi gas at unitarity calculated within the 
model of fermionic and pair excitations with parameters determined using the path-integral GPF 
method (solid curve) and the NSR theory (dashed curve). Thin dotted curves correspond to the 
temperature range above T c . 



and use the grand-canonical thermodynamic potential given by Eq. < f20t with (Q} 
and (l2TT l. 

The second sound velocity U2 characterizes the temperature waves in which 
the motion of the normal and superfluid fractions is out-of -phase. It is determined 
by the formula 



The formulae d26l) and d28t are valid as far as the first-sound and second-sound 
modes are decoupled. Following Ref. — , we assume that the above condition is 
fulfilled for a cold Fermi gas. 

In Fig. 5, the first and second sound velocities (divided by the zero-temperature 
Fermi velocity vf ) are plotted as a function of temperature. They are calculated 
within the model of fermionic and pair excitations using the parameters of pair 
excitations determined by the path-integral GPF and NSR methods. 

The first and second sound velocities obtained using the model of fermionic 
and pair excitations are in a reasonable agreement with the results of the analysis 
based on the full NSR thermodynamic s 31 i 32 . In the zero-temperature limit, the 
first sound velocity tends to the same limit as the sound velocity parameter c for 
pair excitations, which is extremely close to the Monte Carlo dat a ' ' 16 ' 17 ' 18 . 
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Fig. 5 (color online) First and second sound velocities for a unitary Fermi gas calculated within 
the model of fermionic and pair excitations with parameters determined using the path-integral 
GPF method (solid curves) and the NSR theory (dashed curves). The symbols show the Monte 
Carlo results for the first sound velocity. Thin dotted curves correspond to the temperature range 
above T r . 



4 Conclusions 



In the present work, the GPF modification • of the NSR scheme has been for- 
mulated in the path-integral representation and extended to the case of imbalanced 
fermions. Within this path-integral GPF approach, we have analytically derived 
the spectra of low-lying pair excitations of the imbalanced Fermi gas with s-wave 
pairing at finite temperatures and extracted the parameters c and A for the pair 
excitations from these results. Using these spectra, the finite-temperature thermo- 
dynamics of the Fermi gas in the superfluid state has been analyzed. The obtained 
internal energy demonstrates a good agreement with the Monte Carlo results and 
is remarkably close to the experimental data for the Fermi gas at unitarity. The 
zero-temperature value of the first sound velocity is in a good agreement with the 
results of the Monte Carlo simulations. The present method allows us to obtain 
the spectra of the elementary excitations and, consequently, the thermodynamic 
parameters of the state for an arbitrary scattering length, at non-zero temperatures, 
and for non-zero imbalance. 
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